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Abstract 

Systems subjected to holonomic constraints follow quite complicated dynamics that could 
not be described easily with Hamiltonian or Lagrangian dynamics. The influence of 
| holonomic constraints in equations of motions is taken into account by using Lagrange 

multipliers. Finding the value of the Lagrange multipliers allows to compute the forces 
induced by the constraints and therefore, to integrate the equations of motions of the 
system. Computing analytically the Lagrange multipliers for a constrained system may 
be a difficult task that is depending on the complexity of systems. For complex systems, 
it is most of the time impossible to achieve. In computer simulations, some algorithms 
using iterative procedures estimate numerically Lagrange multipliers or constraint forces 
by correcting the unconstrained trajectory. In this work, we provide an analytical com- 
putation of the Lagrange multipliers for a set of linear holonomic constraints with an 
q | arbitrary number of bonds of constant length. In the appendix of the paper, one would 

find explicit formulas for Lagrange multipliers for systems having 1, 2, 3, 4 and 5 bonds 
of constant length, linearly connected. 
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I. Introduction. 



The characteristic time scales associated with intramolecular motions are around 10 to 
100 times shorter than the charecteristic time scales associated with translational and 
rotational degrees of freedom of the molecule. In Molecular Dynamic computations a 
convenient way to handle the multiple-time-scale in equation of motions is to treat cova- 
lent bonds between atoms as rigid, this procedure introduces constraints on equations of 
motions ; the use of the Lagrange multipliers method allows to integrate the equations 
of motion [1]. 

From a technical point of view, the Lagrange multipliers method is used in SHAKE and 
RATTLE algorithms : the constrained trajectories of the systems are computed from 
unconstrained trajectories, by using an iterative algorithm that allows to estimate nu- 
merically the Lagrange multipliers [2, 3, 4]. These algorithms are very efficient and quite 
convenient. Nevertheless, a closed analytical computation of the Lagrange multipliers 
is interesting on its own on a theoretical ground and can also be used to improve the 
iterative estimations done in SHAKE and RATTLE algorithms. 

The equations of motions of the constrained dynamics are derived as follow. For a set of 
holonomic constraints, noted as a a = 0, one defines a constrained Lagrangian L c [q,q] as 

L c [q,q] = L[q,q] - ^ X a a a (q,q) (1) 

a 

where L[q, q] is the Lagrangian of the unconstrained systems and X a the set of Lagrange 
multipliers. The equations of motions are given by 



(2) 



dt dq dq 



which may also be written as 
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■q\ = F t + Y / = ^ + E G ^ ( 3 ) 



where the constraint forces are defined by Gi- a . Since holonomic constraints are con- 
served quantities, by requiring that the second time derivatives of all a a vanish, we 
obtain a set of closed equations that allows to compute the Lagrange multipliers. We 
have 

da'rv 



(4) 

= K a + C a + E Z af3^l3 = 



2 



where we have set 

X 

K a = V — Fi ■ Vi(T a 

% 1 

i,j 

/ — 1 mi 

Then, Lagrange multipliers are given by finding the inverse of the Z-matrix, as 

A /3 = -^(Z- 1 ) /3Q (K Q + C a ) (6) 

a 

Obviously, the structure of the Z-matrix is closely related to the geometrical features of 
molecules, its computation can be quite complicated and is depending on the system. 
In this work, by computing the inverse of the Z-matrix, we provide an analytical result 
for holonomic constraints involved in the dynamic of linear molecules (polymers, etc.) 
or having a linear sequence in its structure (branched polymers, etc). For such systems, 
the holonomic constraints can be read as 

<? a {r a ,r a -i) = (r a - r a -if - a 2 = (7) 

where r a are the coordinates of the atom numbered a in the linear sequence of the 
molecule and a the length of the bond between atoms a and a — 1. N is the number of 
holonomic constraints, these constraints involve TV bonds and N + 1 atoms. 
In the next section, we provide a closed analytical computation of the Z _1 -matrix for 
the linear holonomic constraints given by Eq.(7). 



II. Analytical computation of Lagrange multipli- 
ers for linear holonomic constraints. 

We set r a — r a ^\ = au a , where u a is the unitary bond vector linking the atom numbered 
a — 1 to atom a ; the linear holonomic constraints Eq.(7) gives 

Via a = 2a(5i :a - 5i^-i)u a (8) 

where we use the Kronecker symbol S n:P . The force acting on atom i due to all holonomic 
constraints is Gj = J2a Gi;a and we found easily that 

E G * = ° 
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This may be easily verified for linear holonomic constraints. 

Thus, for a set of linear holonomic constraints, the Z-matrix is given by 

2a? 

Z af3 = — Z^($i,a ~ $i,a-l){&i,P ~ $i,p-l)u a Up 



4a 



(9) 



In the following, we set 

la = ^^a-i • u a and Z aj3 = -i a S a ,p+i + 8 a ,p - jp8 a ,p-i (10) 

In appendix, computations for N = 1,2,3,4 and 5 are given explicitly. With the linear 
holonomic constraints the Z-matrix is a banded symmetric matrix. To compute the 
inverse matrix we take into account the properties of Z. The inverse matrix is built by 
a sequence of multiplications on the right and the left in such a way that the identity 
matrix appears progressively on the diagonal. More precisely, according to Eq.(10), if 
we perform the transform 

Z — ► B\A\ZA X B X = Z^ (11) 
with matrix A\ and B\ defined by 

(Ai)ij = Sij + 72^,i5 2 j and (B^ = S id + ( — = - l)<^2<5 2 ,j (12) 

V 1 - T2 

then, the 2x2 identity matrix appears in for 1 < % < 2 and 1 < j < 2. The 
matrix multiplications by A\ and A\ result in vanishing the non-diagonal coefficient (i.e. 
Zyi = Z>2\ = —72), while matrix multiplications by B\ and B\ scale to 1 the diagonal 
coefficient (A\ZA{)22- The matrix Z^ is similar to the original matrix Z, with 72 and 
73 transformed as 

72 — ► 72 = 

, 73 (13) 

V 1 - T2 2 

because of Eq.(ll). 

We are then in a position to build the identity matrix from Z, by transformations similar 
to Eq.(ll). The matrix Z^ is computed from Z*™ -1 ) by making 

ZW = B^AiZ^AA (14) 

With a convenient choice of the matrix A n and B n , we will finally obtain Z^- 1 ) = I N 
and therefore, we will have built a matrix C such as 

N-l 

C l ZC = I N with C= Yl (A n B n ) (15) 

n=l 
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and thus 

Z~ l = CC l (16) 

To this end, taking into account the structure of Eq.(13) for j' 3 , we define the sequence 
Q n by 

^ = 1-tJ- (17) 

with 71 = 0, ill = 1 and, for convenient reasons, 7^+1 = 0, that gives &n+i = 1- With 
these notations the matrix A n and B n are given by 

(18) 



(B n )ij =S id + (-—-l) 



^,(n+l)<^(n+l),j 



From the definition of matrix A n and B n and with Eq.(16), we may easily compute 
determinants as 

Det(4») = l ; Det(B n ) 



f2 n+ i 
and 

D ^ J UAr^ <19> 

Det(Z) is closely related to the metric determinant that is used in computations of 
ensemble averages [5, 6, 7]. 

Eqs.(16) and (17) allow to compute the Z _1 matrix coefficients, the matrix is symmetric 
Z~ x )ij = Z~ 1 )ji, and after some algebra we found 

e n (jp— 1: 

"i j,=i+im=!+i "i 



(20) 



1 j jV p -1 

fo-<, : = <U n Mi+ e n (7^-1)] 



U fi i n=i+l ^™ p =j+l m =j+l 



With Eqs.(10) and (20), we may verify that Z~ X Z = I N . Eqs.(19) and (20) are the main 
results of this paper. 

According Eq.(6), we have to compute K a and C a to achieve the analytical computation 
of the Lagrange multipliers. With Eqs.(5) and (8), we found 

K a = —(F a -F a _ 1 )-u a (21) 
m 



5 



and 



Col 



;(p a -Pa-iY 



(22) 



with the momentum of atoms, p a = mv a = mr a and notations defined previously. 
Lagrange multipliers are then given by applying Eq.(6). 

So far, we have considered that all atoms, involved in the set of holonomic constraints, 
have the same mass m and the length of each bond is equal to a. We consider now that 
the sequence of atoms mass is {m<j}ie[o,7V] and the length of bonds are {ai}ie[i,7V]- Then, 
the momentum of atoms is now p a = m a v a = m a r a and Eq.(8) becomes 



Therefore, the Z-matrix is now given by 



(23) 



L0 n 



2ai 



m a _i 
1 



(24) 



+ 

m a -i m. 



-!-) 

771 / 



Eq.(24) is similar to Eq.(9). We may transform Eq.(24) to Eq.(lO) with the help of a 
diagonal matrix D, according to 

Z = DZD (25) 

with 



D 



Then, Z is again given by Eq.(lO), but with 

la - 



m a m a -2 



(m a -i + m a ){m a -i + m a „ 2 ) 



(u a _i • u a ) 



(26) 



(27) 



where inertia parameters of bonds appear explicitly (the canonical partition function 
of freely jointed chains is depending on inertia parameters of the bonds as defined in 
ref.[8]). One may note also that if m a = m, Eq.(27) corresponds to Eq.(lO). 
The inverse of the Z-matrix is now given by 



Z- 1 = DZ- l D = D(CC l )D 
with Eq.(27), (17) and (20) ; and K a and C a are given by 

K a = 2a a {— — — -) ■ u a 

and 



•m a m a -i' 



Pa-l 

m a -i 



(28) 

(29) 
(30) 
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III. Discussion. 



From an analytical point of view, an interesting result is given by the computation of 
the determinant of Z (Eq.(19)). For holonomic constraints, this determinant is known 
as the metric determinant ; and for non-hamiltonian dynamical systems [5, 6, 9] it can 
be related to the Jacobian determinant, some ensemble averages can be formulated with 
the help of this determinant [5, 6, 10]. According to computations done in ref.fll], we 
obtain 

i JY l 

P n W^et^" 1 ) = P II &ln^— 

J ° n=2 ^ n=2 fi^Z) ^ 

- [iv-i] He [iv]/i. 3. 1 F 
l 2'2'4'"''4 j 

where ^"^He^ is the multiple hyper geometric function defined in ref.[8], related to the 
canonical partition function of freely jointed chains. In particular, for N = 2, Eq.(31) is 
equivalent to the elementary relation 

f 5 d-io „ ,1 1 3 1 . 1 , , . 

2/ - r 3= = 2 F 1 =2arcsiii- (32) 

Jo A /l- 7 2 2 2 2 4 2 



In ref.[8], we have shown that the multiple hyper geometric function [ Ar_1 lHe^ is defined 
for any mass sequence, Eq.(31) may therefore be extended to any mass sequence, and 
the multiple hyper geometric function ^"^He^ must be evaluated at a point defined 
by the inertia parameters of the coupling between bonds [8]. The inertia parameter of 
the coupling between bonds a and a — 1 is given by 

m a m a - 2 

(m a _i + m a )(m a ^i + m a - 2 ) 

this parameter appears explicitly in Eq.(27), so in the integral too (one may note also 
that if m a = m, for all a, then we have x a = 1/4). For any mass sequence, the relation 
equivalent to Eq.(31) is 

/ dm... / dTivVDet^- 1 ) = ( I] 4 /2 ) [N - 1] ^ N ^p{x i } ie[2 ,N ] ) (34) 
Jo Jo n=2 

with 7 a given by Eq.(27). Eq.(34) can be considered as an integral definition of the 
multiple hypergeometric function ^"^He^ ; this multiple hyper geometric function is 
quite complicated, an iterative approximation scheme, called Independent Motions Ap- 
proximation (IMA), is described in ref.[8]. 

From the point of view of the physical chemistry, the main purpose of this paper was to 
compute analytically the Lagrange multipliers for a set of linear holonomic constraints in 
view of applications to molecular simulations or statistical physics of complex systems. 
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For applications of these analytical results in molecular simulation, the main drawback is 
that the number of coefficients in the Z _1 -matrix, that have to be computed and stored, 
scales as N 2 (the matrix is fully filled, see appendix for some examples), N being the 
number of bonds linearly connected in the molecule. Therefore, for practical implemen- 
tations, this brute force method of computing Lagrange multipliers is certainly not as 
efficient as SHAKE or RATTLE algorithms, at least when N is rather large. Another 
point that one should kept in mind in trying to implement these analytical results in a 
computer code is that, because of C a , velocities (or momentum) contribute to Lagrange 
multipliers. This may rise to some technical problems in the Verlet- velocity algorithm 
and in constant temperature Molecular Dynamics computations, since forces at t + St 
have to be computed to obtain the velocity at t + St [1, 2, 4]. Nevertheless, clever uses 
of these analytical results may certainly be used to improve the efficiency and accuracy 
of RATTLE algorithms ; for instance, when many bonds verify u a -i ■ u a ~ 1. 
On general grounds, despite the complexity of analytical formulas of Z~ 1 given by 
Eq.(20), we believe that the analytical results presented in this work are of some in- 
terest, since matrix similar to Z appears frequently in many models or theories in one 
dimension or with nearest neighbors interactions (spin models, harmonic chains, etc.). 
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Appendix : Lagrange multipliers for N = 1,2,3,4 and 5. 



In this appendix, we give explicitly the analytical formulas for linear holonomic con- 
straints that correspond to N=l, 2, 3, 4 and 5 bonds. To obtain analytical results for 
the matrix Z~ l , we compute firstly from the recursive relation in Eq.(17), 



r n| 



l-7l 



l-7l 



l- 72 2 -73 - 74 2 (l- 7 2 2 ) 
1-71-7 2 

1 - 72 - 7 3 2 - 7l(l - 7 2 2 ) - 7s(l - 7 2 2 - 7 3 2 ) 



l- 7 2 2 - 7 3 2 - 74 2 (l- 72 2 ) 
The computation of Z~ x for N = 1, 2, 3, 4 and 5 is done by using Eq.(20). 

• N = 1 

For N = 1, Z = Z- 1 = (1) and, following Eqs.(20) and (21), we have 



2a 2 o 

-(F 1 - F ) ■ in and C x = —( Pl - p ) 2 



m 



Then, the Lagrange multiplier is 



\ 1 = U( Fl -F 



O) ■ U! + 



1 



(Pi 



2a L ' ma 
the constraint force acting on atom numbered 1 is given by 



-Po) 2 ] 



d = -Ax Vi£7i = - [(Fx - F ) ■ «i + — (pj - P ) 2 ]ui 
and the force acting on atom is Go = — G\. 



ma 



• N = 2 

For N = 2, we have 



1 

"72 



72 



and iT 1 



1 / 1 72 

1 - 72 V 72 1 



(A.l) 



(A.1.1) 



(A.1.2) 



(A.1.3) 



(A.2.1) 
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therefore, with 72 = iii ■ M2/2, Lagrange multipliers are given by 
1 



A, 



2a(l 



(Ml-M 2 )S 
4 / 



[(Fi - F ) • ui + — (Pi-Po) 2 
L ma 

+^(«l • « 2 )[(F 2 - Fi) • «2 + ^;(P 2 - Pi) 2 ]] 



ma 



A 2 = 



2a(l 



(Ml-tX 2 ) 2 



[-(til • u 2 )[(Fi - F ) • ui + — (Pi - p ) 2 ] 



(A.2.2) 



ma 



+ (F 2 -F 1 )-u 2 + —(p 2 - Pl ) 2 ] 



ma 



The forces acting on atoms are 



Go = —2a\\Ui 
G\ = 2aAiiti — 2a\2U2 
G2 = 2a\2U2 



• N = 3 

With iV = 3, matrix Z and Z _1 are given by 



and 



1 -72 
-72 1 -73 
-73 1 



/ 1 - 73 72 7273 \ 
72 1 73 



1 - 72 - 7s 



(A.2.3) 



(A.3.1) 



(A.3.2) 



V 7273 73 1 - 72 / 
As for N = 1, following Eqs. (20) and (21), for 1 < i < 3, K- L and Cj are given by 



2a 2 
Ki = —{Fi- Fj_i) • and d = —{ Pi - p^y 



m 



m^ 



(A.3.3) 
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Then using Eqs.(6), (A. 3. 2) and (A. 3. 3) we may compute Lagrange multipliers and forces 
acting on atoms 



(A.3.4) 



' Go 


= — 2aAiiti 






= 2aX\Ui — 


2aX 2 u 2 


< 

G 2 


= 2a\ 2 U2 - 


2aX 3 ii 3 


> G 3 


= 2aX 3 u 3 





The complicated part of the force being held in Lagrange multipliers. 

• N = 4 

For N = 4, we have 





1 


-72 








\ 




-72 


1 


-73 












-73 


1 


-74 




V 








-74 


1 


/ 



(A.4.1) 



and the matrix coefficients of Z 1 are given by 



z- 1 = 



!-72 -7a -iiiX-iD 



( 1 - 74 - 73 72(1 - 74 ) 7273 727374 \ 

1-74 73 7374 

l-7 2 2 (1 -71)74 

\ 1 - 72 - 7 3 2 J 



Ki and C« are given by Eq.(A.3.3), but 1 < i < 4 and forces by 



' Go 


= — 2aAiUi 






= 2aXiiii — 


2aX 2 U2 


G 2 


= 2aX 2 U2 - 


2aX 3 u 3 


G 3 


= 2aX 3 u 3 - 


2aA4tt4 


, G 4 


= 2aA4tt4 





(A.4.2) 



(A.4.3) 



where Lagrange multipliers are computed by using Eq. (A.4.2). 
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For N = 5, the Z-matrix is 



Z = 



and after some algebra using Eq.(20), we find 



/ 


1 


-72 











\ 




-72 


1 


-73 















-73 


1 


-74 















-74 


1 


-75 




V 











-75 


1 


/ 



(A.5.1) 



1 



1 - 72 - 7s - 74(1 - 72) - 75(1 - 72 - 7s) 



(A.5.2) 



1-75 -74 72(1 - 75 - 7l) 7273(1 - 75) 727374 
- 7 3 2 (l- 7 5 2 ) 

1-75 -74 73(1-71) 7374 



72737475 



737475 



(l- 72 2 )(l- 7 5 2 ) 74(1 - 7 2 2 ) 7475(1 - 7 2 2 : 



1 - 72 - 73 2 75(1 - 72 - 7 3 2 ^ 



2 2 
"72-73 

7l(l- 72 2 ) 



V 



Then, Lagrange multipliers are obtained as previously. 

The computations done for N = 4 and 5 show how the complexity of Lagrange mul- 
tipliers grows with the number of holonomic constraints. On Eqs.(A.4.2) and (A.5.2), 
the equivalence of both directions of labelling atoms, for linear holonomic constraints, 
appears explicitly in the structure of the matrix. 
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